Remote synchronization reveals network symmetries and functional modules 
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We study a Kuramoto model in which the oscillators are associated to the nodes of a complex 
network and the interactions include a phase frustration, thus preventing full synchronization. The 
system organizes into a regime of remote synchronization where pairs of nodes with the same net- 
work symmetry are fully synchronized, despite their distance on the graph. We provide analytical 
arguments to explain this result and we show how the frustration parameter affects the distribution 
of phases. An application to brain networks suggests that anatomical symmetry plays a role in 
neural synchronization by determining correlated functional modules across distant locations. 



PACS numbers: 05.45.Xt, 89.75. Fb, 89.75.Kd 

Synchronization of coupled dynamical units is an ubiq- 
uitous phenomenon in nature [1 . Remarkable examples 
include phase locking in laser arrays, rhythms of flashing 
fireflies, wave propagation in the heart, and also normal 
and abnormal correlations between different regions of 
the human brain [2r[4]. In 1975 Y. Kuramoto proposed a 
simple microscopic model to study collective behaviors in 
large populations of interacting elements [5] . In its origi- 
nal formulation the Kuramoto model describes each unit 
of the system as an oscillator with a pre-assigned natural 
frequency, which interacts with all the other units and 
continuously readjusts its frequency in order to minimize 
the difference between its phase and the phase of all the 
other oscillators. This model has shown very success- 
ful in understanding the emergence of synchronization 
and, over the years, many variations have been consid- 
ered [6HS]. Recently, the Kuramoto model has been also 
extended to the case where the oscillators are the nodes of 
a complex network [2J |H1 [TO] , and it has been found that 
the topology of the interaction network has a fundamen- 
tal role in the emergence and stability of synchronized 
states [11] , In particular, the presence of modular struc- 
tures in the network has a relevant impact on the path to 
synchronization [12-16 , so that practice, units that are 
spatially close to each other on the network, or belong 
to the same tightly connected group of nodes (module 
or community) j!7j , have a higher chance to have similar 
dynamics. This implies that nodes in the same struc- 
tural module often share similar functions, which is a 
belief often supported by empirical findings [3 . How- 
ever, various examples are found in nature where func- 
tional similarity is instead associated with morphological 
symmetry. In these cases, units with similar roles, which 
could potentially swap their position without altering the 
overall functioning of the system, appear in remote lo- 
cations of the network. Some examples include corti- 



cal areas in brains [18], symmetric organs in plants and 
vertebrates, and even atoms in complex molecules |19j . 
Therefore, identifying the sets of symmetric units of a 
complex system might be helpful to understand its orga- 
nization. Finding the global symmetries in a graph, the 
so-called graph automorphism problem, is a well studied 
problem in graph theory, but to date it is still unknown if 
constructing the automorphism group of a graph is poly- 
nomial or NP-complete, even if there exist polynomial- 
time algorithms for graphs with bounded maximum de- 
gree [20] . Recent works have focused instead on the prob- 
lem of defining and detecting local symmetries in complex 
networks [2TJ H2J. Nevertheless, the interplay between 
the structural symmetries of networks and the dynami- 
cal properties of processes occurring over them has been 
studied only marginally [23j . 

In this Letter we show that network symmetries play a 
central role in the synchronization of a system. We con- 
sider networks of Kuramoto oscillators having identical 
natural frequencies, in which a phase frustration param- 
eter forces connected nodes to maintain a finite phase 
difference, thus hindering the attainment of full synchro- 
nization. We prove that the configuration of phases at the 
stationary state depends on the symmetries of the under- 
lying coupling network. In particular, two nodes with the 
same symmetry have identical phases, i.e. are fully syn- 
chronized, and this happens despite the distance of the 
two nodes on the graph. Such a remote synchronization 
behavior is here induced by the network symmetries and 
not by an initial ad-hoc choice of different natural fre- 
quencies [23], and can therefore be used as a method to 
find symmetric nodes in large complex networks. 

In our model N identical oscillators are associated to 
the nodes of a connected graph G(Af,£), with N = \J\f\ 
nodes and K — \C\ links. Each node i is characterized, at 
time t, by a phase 0i(t) whose time evolution is governed 
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by the equation: 
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where uj is the natural frequency, identical for all the os- 
cillators, and A = {a^ } is the adjacency matrix of the 
graph. The model has two control parameters: A > 
accounting for the strength of the interaction, and a, 
the phase frustration parameter ranging in [0, 7r/2[. Ku- 
ramoto systems with frustration and heterogeneous cou- 
pling strengths have been employed to study the emer- 
gence of chimera states, i.e. of stationary states in which 
only a subset of the oscillators attain synchronization 
while the others keep drifting incoherently [Ml US] . 

When a = 0, our model reduces to a network of iden- 
tical Kuramoto oscillators [8] . This is different from 
the classic case of non-identical oscillators, where the fre- 
quency distribution tends to separate the phases and, as 
a result, there is a transition from an incoherent state 



equal to 0) to 



(with order parameter r = 

a synchronized one (r =^ 0) at a critical value A c of the 
coupling, with tightly connected nodes attaining similar 
phases for smaller values of A [13j [Uj . In fact, when all 
the oscillators have the same natural frequency, the only 
attractor of the dynamics is the fully synchronized state, 
no matter the value of A, or the structure of the graph. 
Therefore, starting with random initial phases, the sys- 
tem will reach a fully synchronized state (r — 1) after 
a transient dynamics which reflects the structure of the 
graph, so that nodes belonging to the same structural 
module evolve similarly in time [12] . In conclusion, when 
a = 0, in both the case of identical and non-identical 
oscillators, the nodes with similar phases will be those 
belonging to the same structural module. 

Instead, the introduction of a phase frustration a^0 
forces directly connected oscillators to maintain a con- 
stant phase difference [2B]. In particular, we found that 
for a sufficiently small the dynamics of the model reaches 
a stationary state where the oscillators at two nodes with 
the same symmetry have exactly the same phase, and 
such a phase differs from the phases of nodes with differ- 
ent symmetries. Let us first illustrate this behavior, and 
the effect of the frustration parameter a with the three 
graphs G a , Gb and G c shown in Fig. [I] In Fig. [2] we re- 
port the results of the numerical integration of Eqs. Q 
on the graph G a . In the three topmost panels of Fig. [2] 
we plot the time evolution of the phases of the seven 
nodes for three values of a. We find that, after a tran- 
sient, the system settles into a stationary state in which, 
at any time t, the phases are grouped into four differ- 
ent trajectories: 6»i(t), 9 2 {t) = 9 3 {t), 6» 4 (<) = 7 {t) and 
85(1) = &6(t)- This happens for any value of < a < tt/2, 
with the only difference that by increasing a we better 
separate the four trajectories as shown in the first three 




FIG. 1. (color online) The presence of frustration reveals 
clusters of symmetric nodes. Reported with a color code are 
the phases of the nodes at a given time in the stationary state, 
(a) In the first graph (G a ), node 2 is synchronized to node 
3, node 4 to node 7, and node 5 to node 6. (b) In a finite 
chain (Gb), pairs of nodes symmetrically placed with respect 
to the central node are perfectly synchronized, c) In a finite 
Bethe lattice (G c ) all the nodes placed at the same distance 
from the center have equal phases. 



panels. This is confirmed by panel d), which shows the 
standard deviation of the oscillator phases as a function 
of a. The corresponding four clusters of nodes are identi- 
fied by reporting with a color-code in Fig. [T] the values of 
the phases of the seven nodes at a given time t. We no- 
tice that each cluster groups together all the nodes with 
the same symmetry. In this way two distant nodes of the 
graph, such as for instance node 4 and node 7, are fully 
synchronized even if the other nodes in the paths con- 
necting them have different phases. In this respect, what 
we observe is a remote synchronization |23) induced by 
network symmetries. We have found similar results for 
the linear chain and for the Bethe 3-lattice (see nodes 
with the same colors in Gb and G c in Fig. [fj. 

We can explain these numerical results as follows. If 
a is very small, the difference 9i — 9j is small for any 
couple of nodes i and j, and Eq. ([I]) can be linearized. 
We obtain: 



= d — A 



A' 



(2) 



where are the entries of the Laplacian matrix of the 
graph L = D — A, and D is a diagonal matrix whose 
non-zero elements are the node degrees. Without loss of 
generality, we can set A = 1, Li = , and we can assume 
that at the stationary state 9i = f2,Vi, so that the phases 
must satisfy the equations Sj=i -^ij^j = a [(&) ~ &»] a ^ 
any time, or equivalently: 



L6 = a[(k)l-k] 



(3) 



In a connected graph the Laplacian matrix has one null 
eigenvalue and the system of Eqs. ([3| is singular, i.e. 
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FIG. 2. (color online) After an initial transient the system 
is phase-locked in a stationary state where nodes with the 
same symmetry have the same phase. The figure refers to 
the graph G a in Fig. The three top panels report three 
different values of the frustration parameter, respectively a) 
a = O.f, b) a = 0.5, c) a = 0.8, while in panel d) we show 
how the spreading of the phases ag increases with the value 
of Q. 



it has one free variable. Consequently, at each time t 
we can solve the system by computing the phase dif- 
ference between each node and a given node chosen as 
reference. For instance, for G a , we can define the vari- 
ables <f>j(t) = 9j(t) — Oi(t),j — 2,..., 7, and by solv- 
ing the equations we obtain: 4>2 = 4>z = a l(k) — 2], 
04 = 7 = 2a[(fc) - 2], and 5 = </> 6 = 3a[{k}-2]. 
This is in agreement with the results of the simulations: 
in G a the phases are clustered into four groups, with 
nodes with the same symmetry having the same phase, 
and nodes with different symmetries being separated by 
a phase lag that depends on a as in the relations found 
above. We can derive analogous analytical expressions 
for chain graphs and Bethe lattices. For instance, for 
graph Gb in Fig[T]we obtain: 



-n-H 



(4) 



If we solve Eq. Q for i — 77-, we obtain 9 n — 0, 



^{n+l) 
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a. Consequently, two nodes 



symmetrically placed with respect to node will have 
identical phases. 

We now provide a general argument to explain why 
the synchronization in the frustrated Kuramoto model is 
related to the graph symmetry. A graph G(Af, C) has a 
symmetry if and only if it is possible to find a bijection 
7r : Af — > Af which preserves the adjacency relation of 
G, i.e. which is an automorphism for G. Formally, this 
means that there exists a permutation matrix P = P(tt), 
i.e. a matrix with each row and each column having ex- 
actly one entry equal to one and all others equal to 0), 



such that PAP^ 1 = A. Since P is a permutation matrix, 
the matrix product PA swaps pairs of rows of A, while 
AP^ 1 swaps pairs of columns of A. Therefore, if P cor- 
responds to an automorphism of G then PAP -1 actually 
performs a relabeling of the nodes of the original graph 
which preserves the adjacency matrix. This implies that 
P commutes with A, i.e. PA = AP. In general a graph 
can admit more than one automorphism. For instance, 
graph G a in Fig|T]has at least three non-trivial bijections 
which preserve the adjacency matrix, namely: 



7Tl 
7T3 



(1,2,3,4,5,6,7) 
(1,2,3,4,5,6,7) 
(1,2,3,4,5,6,7) 



-> (1,3,2,4,5,6,7) 
(1,2,3,7,6,5,4) 
-> (1,3,2,7,6,5,4) 



Node 2 and node 3 are symmetric because it is possible 
to find a relabeling of the nodes of G a (e.g. by means 
of either m or Tts) so that node 2 is mapped into node 
3 and vice-versa, and the adjacency matrix of G a is left 
unchanged. Similarly, for the pairs {4,7} and {5,6} we 
have two different relabelings (i.e. Tt2 and ^3) mapping 
one node of each pair on the other and preserving adja- 
cency relations. Therefore, in terms of symmetries, the 
graph G has four different classes of nodes: C\ = {1}, 
C 2 = {2,3}, C 3 = {4,7}, C 4 = {5,6}. Now, if a per- 
mutation of the nodes of a graph G is an automorphism, 
then PLP- 1 = PDP- 1 - PAP- 1 =D-A = L, i.e. the 
associated permutation matrix P also commutes with the 
Laplacian matrix of the graph. If we now left-multiply 
both sides of Eq. §byPwe get PL0 = aP [(k)l - k]. 
Since P commutes with L and permutes symmetric nodes 
(notice that symmetric nodes have the same degree) , then 
PL = LP and Pk = fc, so we obtain: 



LP0 = a[{k)l-k] 



(5) 



Making use of Eq. ([3| and Eq. ([5| we get the linear sys- 
tem: LP9 = LB. This system is singular, i.e. it has one 
free variable. Again, it can be solved by leaving free one 
of the variables 0j, setting <j>j = 6j — 9i and considering 
the new system LP<p = Lcf>. The matrix P is the matrix 
obtained from P by removing the row and the column 
corresponding to node i. HP does not permute node i 
with another node, then P is still a permutation matrix. 
Analogously, the matrix L is the reduced Laplacian, i.e. 
the non-singular matrix obtained from the Laplacian by 
deleting the i— th row and the i— th column. By left- 
multiplying by i" 1 , which is not singular, we obtain: 



P0 = 4> 



(6) 



Since P(f> is a permutation of the phases of symmetric 
nodes, Eq.(|6]) implies that the phases of symmetric nodes 
will be equal at any time, whereas by solving Eq. ([5| we 
can get the actual values of the corresponding phases. 

Application to the brain. - To investigate the functional 
role of symmetry in the human brain we have consid- 



4 



(a) 





200 t (sec) 



(c) 




log(A0) 



FIG. 3. (color online) (a) Brain sites with similar and dissimilar phases of the frustrated Kuramoto model are colored and 
superimposed onto an anatomical image, (b) Examples of functional data from one subject recorded at the brain areas indicated 
in plot (a). Colors are the same as those used in the anatomical image, (c) Functional correlation Z between pairs of nodes 
as a function of their phase differences A8. Black solid curve corresponds to the average value over all the subjects, while the 
gray area covers the 5 th and the 95 th percentiles of the distribution. The dashed horizontal line indicates the threshold for 
statistical significant correlations (p < 0.05, corrected for multiple comparisons). 



ered anatomical and functional brain connectivity ma- 
trices defined on the same set of iV = 90 nodes repre- 
senting different areas of the brain (see details in Ap- 
pendix). Wc have first constructed a graph where the 
links represent anatomical connectivity as obtained from 
DW-MRI data [57] and used this graph as a backbone 
network to integrate Eq. (JlJ. We identified candidate 
pairs of anatomically symmetric brain areas by means of 
an agglomerative clustering algorithm which groups to- 
gether nodes having close phases at the stationary state 
(full dendrogram is reported in Fig. [4| . Then, in order 
to assess the role of symmetry in the real activity of the 
brain, we have looked at the correlations between the 
BOLD fMRI time-series recorded during resting state on 
healthy individuals. We illustrate the results in Fig. [3] 
Consider nodes 57 and 74 corresponding respectively to 
green and blue areas in panel (a). Not only the two areas 
are spatially separated, but there is no edge connecting 
the two corresponding nodes in the anatomical connec- 
tivity network. However, the two nodes are detected as 
a candidate symmetric pair since at the stationary state 
of the Kuramoto dynamics in Eq. ([!]) the oscillators as- 
sociated to these two nodes have very close phases (see 
dendrogram in Fig. |4|. As shown in Figure [3Jd, we found 
that also the BOLD fMRI signals corresponding to nodes 
57 and 74 are strongly synchronized. We obtain remark- 
ably different results when we consider another pair of 
nodes, such as node 74 and node 76. These two nodes 
correspond to two spatially adjacent areas of the brain 
(the red and blue regions in Fig. [3^,) and are directly con- 
nected in the anatomical connectivity network. However, 
at the stationary state of Eq. (JlJ the phase difference of 
the oscillators associated to node 74 and 76 is quite large. 
Interestingly, in this case the fMRI time-series associated 
to these nodes are much less similar to each other (see 
the two bottom trajectories reported in Fig. pjj: 



To quantify this effect, we have computed the similar- 
ity between the fMRI time scries associated to the nodes 
by means of linear correlation, followed by the Fisher's 
Z-transform (to stabilize the variance of the estimator) 
and a correction for correlated samples. In Fig. [3J: we 
plot the average functional correlation Z between pairs 
of brain areas as a function the phase differences A9 be- 
tween the phases of the corresponding oscillators as ob- 
tained from the dynamics of Eq. ([I]) on the anatomical 
connectivity network. The fact that Z decreases with 
A9 reveals that structural symmetry plays an important 
role in determining human brain functions. In fact, the 
functional activities of anatomically symmetric areas can 
be strongly correlated, even if the areas are distant in 
space. In this way the study of anatomical symmetries 
in neural systems might provide meaningful insights into 
the functional organization of distant brain neural assem- 
blies during diverse cognitive or pathological states [18] . 
Applied to other connectivity networks as a method to 
detect network symmetries, our study could provide new 
insights on the interplay between structure and dynamics 
in complex systems. 
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APPENDIX 

Data acquisition and pre-processing. - The anatomical 
connectivity network is based on the connectivity ma- 
trix obtained by Diffusion Magnetic Resonance Imaging 
(DW-MRI) data from 20 healthy participants, as de- 
scribed in [28] . The elements of this matrix represent 
the probabilities of connection between the 90 anatom- 
ical regions of interest (N = 90 nodes in the network) 
of the Tzourio-Mazoyer brain atlas [29]. These proba- 
bilities are proportional to the density of fibers between 
different areas, so each element of the matrix represents 
an approximation of the connection strength between the 
corresponding pair of brain regions. 

The functional brain connectivity was extracted from 
BOLD fMRI resting state recordings obtained as de- 
scribed in [5U] • All fMRI datasets (segments of 5 minutes 
recorded from 15 healthy subjects) were co-registered to 
the anatomical dataset and normalized to the standard 
MNI (Montreal Neurological Institute) template image, 
to allow comparisons between subjects. As DW-MRI 
data, normalized and corrected functional scans were 
subsampled to the anatomical labeled template of the 
human brain [29]. Regional time series were estimated 
for each individual by averaging the fMRI time series 
over all voxels in each region (data were not spatially 
smoothed before regional parcellation) . To eliminate low 
frequency noise (e.g. slow scanner drifts) and higher fre- 
quency artifacts from cardiac and respiratory oscillations, 
time-series were digitally filtered with a finite impulse re- 
sponse (FIR) filter with zero-phase distortion (bandwidth 
0.01-0.1 Hz) as in [3D]. 

Functional synchrony. - A functional link between two 
time series Xi (t) and Xj (t) (normalized to zero mean and 
unit variance) was defined by means of the linear cross- 
correlation coefficient computed as = (Xi(f)Xj(t)), 
where (•) denotes the temporal average. For the sake 
of simplicity, we only considered here correlations at lag 
zero. To determine the probability that correlation val- 
ues are significantly higher than what is expected from 
independent time series, 7V/(0) values (denoted r^) were 
firstly transformed by the Fisher's Z transform 



rected by the following approximation: 
1 1 2 



(7) 



Under the hypothesis of independence, Zij has a nor- 
mal distribution with expected value and variance 
l/(df — 3), where df is the effective number of degrees 
of freedom [5H433] . If the time series consist of indepen- 
dent measurements, df simply equals the sample size, N. 
Nevertheless, autocorrelated time series do not meet the 
assumption of independence required by the standard sig- 
nificance test, yielding a greater Type I error [31-55]. In 
presence of auto-correlated time series df must be cor- 



df ~ N N 



Xl r «( T ) r Ji( r )' 



(8) 



where r xx (r) is the autocorrelation of signal x at lag r. 

To estimate a threshold for statistically significative 
correlations, a correction for multiple testing was used. 
The False Discovery Rate (FDR) method was applied 
to each matrix of values [53]. With this approach, 
the threshold of significance Z t h was set such that the 
expected fraction of false positives is restricted to q < 
0.05. 

Clustering of phase values.- To identify brain areas 
related by a topological symmetry, we use the anatomi- 
cal connectivity obtained from the DW-MRI data as the 
connectivity matrix (N — 90 nodes) in Eq. (1). A stan- 
dard hierarchical agglomerative clustering algorithm is 
then used to identify nodes with similar phases [35] . The 
resulting dendrogram is depicted in Fig. [4] 
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